\section{Mathematical Model of a Step Discontinuity}
A step in a radial waveguide can be analyzed by field matching the modes form one waveguide to the modes in the adjacent waveguide.  The basic steps to solving the problem are the following:
\begin{enumerate}
	\item Expand the fields in each region in terms of its Eigen modes
	\item Apply boundary conditions across the step discontinuity
	\item Test with inner product operators and weighting functions
	\item Solve for the scattering matrix
\end{enumerate}
Each of these steps will be explained in more detail.  

\subsection{Expansion of Fields}
The transverse fields in region (1) are expanded as follows.
\begin{align}
\vec{E}_\bot^{[1]}(\rho,\phi,z) &= \sum_m \vec{e}_m^{[1]}(\phi,z)v_m^{[1]}(\rho)\label{eqn:rwexpE1}\\
\vec{H}_\bot^{[1]}(\rho,\phi,z) &= \sum_m \vec{h}_m^{[1]}(\phi,z)i_m^{[1]}(\rho)\label{eqn:rwexpH1}
\end{align}
Similarly in region (2) they are expanded as,
\begin{align}
\vec{E}_\bot^{[2]}(\rho,\phi,z) &= \sum_n \vec{e}_n^{[2]}(\phi,z)v_n^{[2]}(\rho)\label{eqn:rwexpE2}\\
\vec{H}_\bot^{[2]}(\rho,\phi,z) &= \sum_n \vec{h}_n^{[2]}(\phi,z)i_n^{[2]}(\rho)\label{eqn:rwexpH2}
\end{align}
Note that the number of modes in region 1 does not have to equal the number of modes in region 2, however the number of modes in the larger guide should be greater than the number of modes in the smaller guide in order for the solution to converge \cite{Eleftheriades:1994}.  

\subsection{Application of Boundary Conditions}
The boundary conditions across the step from region 1 to region 2 are,
\begin{align}
\hat{\rho}\times\vec{E}_\bot^{[1]} = \hat{\rho}\times\vec{E}_\bot^{[2]}\label{eqn:rwbce}\\
\hat{\rho}\times\vec{H}_\bot^{[1]} = \hat{\rho}\times\vec{H}_\bot^{[2]}\label{eqn:rwbch}
\end{align}
Substituting (\ref{eqn:rwexpE1}) and (\ref{eqn:rwexpE2}) into (\ref{eqn:rwbce}) and substituting (\ref{eqn:rwexpH1}) and (\ref{eqn:rwexpH2}) into (\ref{eqn:rwbch}) yields,
\begin{align}
\hat{\rho}\times\sum_m \vec{e}_m^{[1]}(\phi,z)v_m^{[1]}(\rho) = \hat{\rho}\times\sum_n \vec{e}_n^{[2]}(\phi,z)v_n^{[2]}(\rho)\label{eqn:rwbce2}\\
\hat{\rho}\times\sum_m \vec{h}_m^{[1]}(\phi,z)i_m^{[1]}(\rho) = \hat{\rho}\times\sum_n \vec{h}_n^{[2]}(\phi,z)i_n^{[2]}(\rho)\label{eqn:rwbch2}
\end{align}

\subsection{Testing with weighting function and inner product}
An inner product for testing or discretizing the above equations is given as,
\begin{align}
\int_{a_1}^{b_1}\int_0^{2\pi}\{\bullet\}\cdot\vec{h}_p^{[1]}(\phi,z)(\rho d\phi)dz\label{eqn:innerprodr1}\\
\int_{a_2}^{b_2}\int_0^{2\pi}\{\bullet\}\cdot\vec{e}_q^{[2]}(\phi,z)(\rho d\phi)dz\label{eqn:innerprodr2}
\end{align}
where the eigen modes from the larger guide is used to match the electric fields and the modes form the smaller guide is used to match the magnetic fields \cite{Eleftheriades:1994}, \cite{Morini:2001}.  This implies that region (1) should be the larger of the two waveguides. Applying (\ref{eqn:innerprodr1}) to both sides of (\ref{eqn:rwbce2})gives,
\begin{multline}
\int_{a_1}^{b_1}\int_0^{2\pi}\hat{\rho}\times\sum_m \vec{e}_m^{[1]}(\phi,z)v_m^{[1]}(\rho)\cdot\vec{h}_p^{[1]}(\phi,z)(\rho d\phi)dz=\\ \int_{a_1}^{b_1}\int_0^{2\pi}\hat{\rho}\times\sum_n \vec{e}_n^{[2]}(\phi,z)v_n^{[2]}(\rho)\cdot\vec{h}_p^{[1]}(\phi,z)(\rho d\phi)dz\label{eqn:rwbce3}
\end{multline}
and applying (\ref{eqn:innerprodr2}) to both side of (\ref{eqn:rwbch2}) yields,
\begin{multline}
\int_{a_2}^{b_2}\int_0^{2\pi}\hat{\rho}\times\sum_m \vec{h}_m^{[1]}(\phi,z)i_m^{[1]}(\rho)\cdot\vec{e}_q^{[2]}(\phi,z)(\rho d\phi)dz=\\ \int_{a_2}^{b_2}\int_0^{2\pi}\hat{\rho}\times\sum_n\vec{h}_n^{[2]}(\phi,z)i_n^{[2]}(\rho)\cdot\vec{e}_q^{[2]}(\phi,z)(\rho d\phi)dz\label{eqn:rwbch3}
\end{multline}
Equation (\ref{eqn:rwbce3}) and (\ref{eqn:rwbch3}) can be manipulated into the following form,
\begin{align}
\sum_mR^{[1]}_{p,m}v_m^{[1]}(\rho) = \sum_nP_{p,n}v_n^{[2]}(\rho)\label{eqn:rwbce4}\\
\sum_mQ_{q,m}i_m^{[1]}(\rho)=\sum_nR^{[2]}_{q,n}i_n^{[2]}(\rho)\label{eqn:rwbch4}
\end{align}
where,
\begin{align}
R^{(1)}_{p,m} &= \int_{a_1}^{b_1}\int_0^{2\pi}\vec{e}_m^{[1]}(\phi,z)\times\vec{h}_p^{[1]}(\phi,z)\cdot\hat{\rho}\;d\phi dz\label{eqn:defR1}\\
R^{(2)}_{q,n} & = \int_{a_2}^{b_2}\int_0^{2\pi}\vec{e}_q^{[2]}(\phi,z)\times\vec{h}_n^{[2]}(\phi,z)\cdot\hat{\rho}\;d\phi dz\label{eqn:defR2}\\
P_{p,n} & = \int_{a_2}^{b_2}\int_0^{2\pi}\vec{e}_n^{[2]}(\phi,z)\times\vec{h}_p^{[1]}(\phi,z)\cdot\hat{\rho}\;d\phi dz\label{eqn:defP}\\
Q_{q,m} &= \int_{a_2}^{b_2}\int_0^{2\pi}\vec{e}_q^{[2]}(\phi,z)\times\vec{h}_m^{[1]}(\phi,z)\cdot\hat{\rho}\;d\phi dz\label{eqn:defQ}
\end{align}
Note that the limits of integration collapsed to the common region over both waveguides.  Also note that the $\rho$ dependence canceled out on both sides and the integrals are only a function of the transverse variables $\phi$ and $z$. 

Equations (\ref{eqn:rwbce4}) and (\ref{eqn:rwbch4}) can be cast into matrix form as,
\begin{align}
\textbf{R}_1v_1 &= \textbf{P}v_2\label{eqn:ecylmatch}\\
\textbf{Q}i_1 &= \textbf{R}_{2}i_2\label{eqn:hcylmatch}
\end{align}
where $\textbf{R}_1$, $\textbf{R}_2$, $\textbf{P}$ and  $\textbf{Q}$ are matrices and $v_1$, $v_2$, $i_1$ and $i_2$ are column vectors.

\subsection{Non-Unifrom Transition (ABCD) Matrix}
\begin{align}
\left[ \begin{array}{c}
v_1 \\
i_1
\end{array}\right] = 
\left[ \begin{array}{cc}
\textbf{T}_{1,1} & \textbf{T}_{1,2} \\
\textbf{T}_{2,1} & \textbf{T}_{2,2} \\
\end{array} \right]
\left[ \begin{array}{c}
v_2 \\
i_2
\end{array}\right]
\end{align}
where,
\begin{align}
\textbf{T}_{1,1} &= \textbf{R}_{1}^{-1} \textbf{P}\\
\textbf{T}_{1,2} &= 0\\
\textbf{T}_{2,1} &= 0\\
\textbf{T}_{2,2} &= \textbf{Q}^{-1} \textbf{R}_{2}
\end{align}


\subsection{Non-Uniform Waveguide Mode Scattering}

Since a scattering matrix is desired, the normalized variables $v(\rho)$ and $i(\rho)$ are expanded as the summation of an incident wave coefficient ($a$) and a reflected wave coefficient ($b$) along with wave propagation functions that describe how and in which direction the wave propagates for a given geometry and coordinate system.  For region 1,  
\begin{align}
v_1(\rho) &= \textbf{v}_1^+(\rho)\;a_1+\textbf{v}_1^-(\rho)\;b_1\label{eqn:v1rho}\\
i_1(\rho) &= \textbf{i}_1^+(\rho)\;a_1+\textbf{i}_1^-(\rho)\;b_1\label{eqn:i1rho}
\end{align}
where for region 2,
\begin{align}
v_2(\rho) &= \textbf{v}_2^+(\rho)\;b_2+\textbf{v}_2^-(\rho)\;a_2\label{eqn:v2rho}\\
i_2(\rho) &= \textbf{i}_2^+(\rho)\;b_2+\textbf{i}_2^-(\rho)\;a_2\label{eqn:i2rho}
\end{align}
where,
\begin{align}
\textbf{v}_k^{\pm} = \textbf{U}_k{v}_k^{\pm}\label{eqn:vectomatvk}\\
\textbf{i}_k^{\pm} = \textbf{U}_k{i}_k^{\pm}\label{eqn:vectomatik}
\end{align}
and where $k=1,2$ and $\textbf{U}_k$ is an identity matrix. A note on convention is that a function $f^+(\rho)$ indicates a wave travelling in the $+\rho$ direction while a $f^-(\rho)$ indicates a wave travelling in the $-\rho$ direction.
Since the geometry of interest is best described in the cylindrical coordinate system and the waves are propagating in the $\hat{\rho}$ direction, then the functions that best describe the wave propagation are Hankel functions and the derivatives of Hankel functions.  The exact definitions are dependent on whether the modes are $TE_z$ or $TM_z$ so they will be kept in general terms for now to derive the scattering matrix. Substituting (\ref{eqn:v1rho})and (\ref{eqn:v2rho}) into (\ref{eqn:ecylmatch}) and substituting (\ref{eqn:i1rho}) and (\ref{eqn:i2rho}) into (\ref{eqn:hcylmatch}) yields,
\begin{align}
\textbf{R}_1\textbf{v}_1^+a_1+\textbf{R}_1\textbf{v}_1^-b_1 &= \textbf{P}\textbf{v}_2^-a_2+\textbf{P}\textbf{v}_2^+b_2\label{eqn:ecylmatch3}\\
\textbf{Q}\textbf{i}^+_1a_1+\textbf{Q}\textbf{i}_1^-b_1 &= \textbf{R}_{2}\textbf{i}_2^-a_2+\textbf{R}_{2}\textbf{i}_2^+b_2\label{eqn:hcylmatch3}
\end{align}
Solving (\ref{eqn:hcylmatch3}) for $b_2$ gives,
\begin{align}
b_2 &= (\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{Q}\textbf{i}^+_1a_1 + (\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{Q}\textbf{i}_1^-b_1 -(\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{R}_{2}\textbf{i}_2^-a_2\label{eqn:hcylmatch4}
\end{align}
Substituting (\ref{eqn:hcylmatch4}) into (\ref{eqn:ecylmatch3}) and rearranging yields,
\begin{multline}
\left[\textbf{R}_2\textbf{v}_1^- - \textbf{P}\textbf{v}_2^+(\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{Q}\textbf{i}_1^-\right]b_1 =\\ [\textbf{P}\textbf{v}_2^+(\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{Q}\textbf{i}^+_1 -\textbf{R}_1\textbf{v}_1^+]a_1 + \textbf{P}[\textbf{v}_2^- - \textbf{v}_2^+(\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{R}_{2}\textbf{i}_2^-]a_2\label{eqn:ecylmatch4}
\end{multline}
If $\textbf{R}_1$ and $\textbf{R}_2$ are identity matrices, in which case they will be if reaction (instead of power) normalization is used then,
\begin{multline}
\left[\textbf{v}_1^- - \textbf{P}\textbf{v}_2^+(\textbf{i}_2^+)^{-1}\textbf{Q}\textbf{i}_1^-\right]b_1 =\\ [\textbf{P}\textbf{v}_2^+(\textbf{i}_2^+)^{-1}\textbf{Q}\textbf{i}^+_1 -\textbf{v}_1^+]a_1 + \textbf{P}[\textbf{v}_2^- - \textbf{v}_2^+(\textbf{i}_2^+)^{-1}\textbf{i}_2^-]a_2\label{eqn:ecylmatch5}
\end{multline}
Solving (\ref{eqn:ecylmatch3}) for $b_1$ gives,
\begin{align}
b_1 &= (\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{P}\textbf{v}_2^-a_2 + (\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{P}\textbf{v}_2^+b_2 - (\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{R}_1\textbf{v}_1^+a_1\label{eqn:ecylmatch6}
\end{align}
Substituting (\ref{eqn:ecylmatch6}) into (\ref{eqn:hcylmatch3}) and rearranging gives,
\begin{multline}
[\textbf{R}_{2}\textbf{i}_2^+-\textbf{Q}\textbf{i}_1^-(\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{P}\textbf{v}_2^+]b_2 =\\ \textbf{Q}[\textbf{i}^+_1-\textbf{i}_1^-(\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{R}_1\textbf{v}_1^+]a_1 + [\textbf{Q}\textbf{i}_1^-(\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{P}\textbf{v}_2^- - \textbf{R}_{2}\textbf{i}_2^-]a_2\label{eqn:hcylmatch5}
\end{multline}
and if $\textbf{R}_1$ and $\textbf{R}_2$ are identity matrices then,
\begin{multline}
[\textbf{i}_2^+-\textbf{Q}\textbf{i}_1^-(\textbf{v}_1^-)^{-1}\textbf{P}\textbf{v}_2^+]b_2 =\\ \textbf{Q}[\textbf{i}^+_1-\textbf{i}_1^-(\textbf{v}_1^-)^{-1}\textbf{v}_1^+]a_1 + [\textbf{Q}\textbf{i}_1^-(\textbf{v}_1^-)^{-1}\textbf{P}\textbf{v}_2^- - \textbf{i}_2^-]a_2\label{eqn:hcylmatch6}
\end{multline}
A generalized scattering matrix can now be defined as,
\begin{align}
\left[ \begin{array}{c}
b_1 \\
b_2
\end{array}\right] = 
\left[ \begin{array}{cc}
\textbf{S}_{1,1} & \textbf{S}_{1,2} \\
\textbf{S}_{2,1} & \textbf{S}_{2,2} \\
\end{array} \right]
\left[ \begin{array}{c}
a_1 \\
a_2
\end{array}\right]
\end{align}
where $b_1$, $b_2$, $a_1$ and $a_2$ are column vectors and $\textbf{S}_{1,1}$, $\textbf{S}_{2,1}$, $\textbf{S}_{1,2}$ and $\textbf{S}_{2,2}$ are partition matrices which are given as,
\begin{align}
\textbf{S}_{1,1} &= \textbf{V}^{-1}[\textbf{P}\textbf{v}_2^+(\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{Q}\textbf{i}^+_1 -\textbf{R}_1\textbf{v}_1^+]\label{eqn:rwgS11}\\
\textbf{S}_{1,2} &= \textbf{V}^{-1}\textbf{P}[\textbf{v}_2^- - \textbf{v}_2^+(\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{R}_{2}\textbf{i}_2^-]\label{eqn:rwgS12}\\
\textbf{S}_{2,1} &= \textbf{W}^{-1}\textbf{Q}[\textbf{i}^+_1-\textbf{i}_1^-(\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{R}_1\textbf{v}_1^+]\label{eqn:rwgS21}\\
\textbf{S}_{2,2} &= \textbf{W}^{-1}[\textbf{Q}\textbf{i}_1^-(\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{P}\textbf{v}_2^- - \textbf{R}_{2}\textbf{i}_2^-]\label{eqn:rwgS22}
\end{align}
where,
\begin{align}
	\textbf{V} &= \left[\textbf{R}_2\textbf{v}_1^- - \textbf{P}\textbf{v}_2^+(\textbf{R}_{2}\textbf{i}_2^+)^{-1}\textbf{Q}\textbf{i}_1^-\right]\label{eqn:auxV}\\
	\textbf{W} &= \left[\textbf{R}_2\textbf{i}_2^+-\textbf{Q}\textbf{i}_1^-(\textbf{R}_1\textbf{v}_1^-)^{-1}\textbf{P}\textbf{v}_2^+\right]\label{eqn:auxW}
\end{align}

\subsection{Propagation Matrices for Radial Waveguides}
For region 1 (\ref{eqn:v1rho}) and (\ref{eqn:i1rho}) can be written in matrix form as,
\begin{align}
\left[ \begin{array}{c}
v_1(\rho) \\
i_1(\rho)
\end{array}\right] = 
\left[ \begin{array}{cc}
v_1^+(\rho) & v_1^-(\rho) \\
i_1^+(\rho) & i_1^-(\rho)
\end{array} \right]
\left[ \begin{array}{c}
a_1 \\
b_1
\end{array}\right]\label{eqn:propmatrix1}
\end{align}
where similarly for region 2,
\begin{align}
\left[ \begin{array}{c}
v_2(\rho) \\
i_2(\rho)
\end{array}\right] = 
\left[ \begin{array}{cc}
v_2^+(\rho) & v_2^-(\rho) \\
i_2^+(\rho) & i_2^-(\rho)
\end{array} \right]
\left[ \begin{array}{c}
b_2 \\
a_2
\end{array}\right]\label{eqn:propmatrix2}
\end{align}

The functional matrix is a function of the independent variable in the direction of propagation.  This matrix will be called the \emph{wave propagation matrix} or just the \emph{propagation matrix} since it describes how the wave propagates for a given coordinate system.  For the case of a radial waveguide, Hankel functions will be used since the the cylindrical coordinate system has been chosen to describe the geometry.  

\subsubsection{Radial Waveguide Propagation Matrix Definitions}
The definitions of the propagation matrices change depending whether the modes are $TE_a$ or $TM_z$.  For $TE_z$ modes the propagation matrix is found to be,
\begin{align}
\left[ \begin{array}{cc}
v_{f}^+(\rho) & v_{f}^-(\rho) \\
i_{f}^+(\rho) & i_{f}^-(\rho)
\end{array} \right] &=
\left[ \begin{array}{cc}
\frac{1}{-k_{f\rho}}{H_{m}^{(2)}}'(k_{f\rho}\rho) & \frac{1}{-k_{f\rho}}{H_{m}^{(1)}}'(k_{f\rho}\rho)\\
{H_{m}^{(2)}}(k_{f\rho}\rho) & {H_{m}^{(1)}}(k_{f\rho}\rho)
\end{array} \right]
\end{align}
For $TM_z$ modes it is found to be,
\begin{align}
\left[ \begin{array}{cc}
v_{a}^+(\rho) & v_{a}^-(\rho) \\
i_{a}^+(\rho) & i_{a}^-(\rho)
\end{array} \right] &=
\left[ \begin{array}{cc}
{H_m^{(2)}}(k_{a\rho n}\rho) & {H_m^{(1)}}(k_{a\rho n}\rho)\\
\frac{1}{-k_{a\rho}}{H_m^{(2)}}'(k_{a\rho n}\rho) & \frac{1}{-k_{a\rho}}{H_m^{(1)}}'(k_{a\rho n}\rho)
\end{array} \right]
\end{align}

\subsubsection{Shift in Reference Plane}
From (\ref{eqn:propmatrix1}) and (\ref{eqn:propmatrix2}) it can be seen that for any given region $k$, $a_k$ and $b_k$ are constants and do not change for the entire region.  This means that once these constants are known, a shift in the reference plane of the scattering matrix can be found by the use of the propagation matrix (\ref{eqn:propmatrix}).  All that needs to be done is that for a new value of $\rho$, a new propagation matrix is calculated and then substituted into the scattering matrix formulas in (\ref{eqn:rwgS11})--(\ref{eqn:rwgS22}).